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It  is  of  great  interest  to  know  if  an  eventual  quantum 
computer  could  solve  a  broad  range  of  hard  “optimiza¬ 
tion”  problems  more  efficiently  than  a  classical  computer. 
An  important  class  is  the  NP-hard  category  [1] 
(nondeterministic-polynomial-time-hard),  for  which  it  is 
believed  that  all  classical  algorithms  take  a  time  which 
grows  exponentially  with  the  problem  size  N. 

The  most  promising  approach  to  solving  optimization 
problems  on  a  quantum  computer  seems  to  be  the  quantum 
adiabatic  algorithm  (QAA),  proposed  by  Farhi  et  al.  [2]. 
The  idea,  which  is  related  to  quantum  annealing  [3],  is  that 
one  adds  to  a  “problem”  Hamiltonian,  3~C P,  whose  ground 
state  represents  a  solution  of  a  classical  optimization  prob¬ 
lem,  a  noncommuting  “driver”  Hamiltonian,  3~C D,  so  the 
total  Hamiltonian  is 

3~C  (s)  =  (1  —  s)J~CD  +  sfhCp,  (1) 

where  .v(/)  is  a  time-dependent  control  parameter.  3~C  P  is 
expressed  in  terms  of  classical  Ising  spins  taking  values 
±1,  or  equivalently  in  terms  of  the  z  components  of  the 
Pauli  matrices  for  each  spin,  fr].  The  simplest  driver 
Hamiltonian  is  then  3~C D  =  —  YJi=  i  0"?  where  erf  is  the 
.r-component  Pauli  matrix. 

The  control  parameter  s(t )  is  0  at  t  =  0,  so  3~C  =  J~C D, 
which  has  a  trivial  ground  state  in  which  all  2N  basis  states 
(in  the  fr  basis)  have  equal  amplitude.  It  then  increases 
with  t ,  reaching  1  at  t  =  'T ,  where  'T  is  the  runtime  of  the 
algorithm,  at  which  point  3~C  =  3~CP.  If  the  time  evolution 
of  s(t )  is  sufficiently  slow,  the  process  will  be  adiabatic. 
Hence,  starting  the  system  in  the  ground  state  of  J~CD,  the 
system  will  end  up  in  the  ground  state  of  3~C P  and  the 
problem  is  solved.  The  time  'T  required  to  find  the  ground 
state  with  significant  probability  is  called  the  complexity. 
The  bottleneck  of  the  QAA  is  likely  to  be  at  one  or  more 
points  where  the  energy  gap  from  the  ground  state  to  the 


first  excited  state  becomes  very  small,  possibly  due  to  a 
quantum  phase  transition. 

Early  numerical  work  [2,4]  on  very  small  systems,  N  < 
24  (for  a  particular  constraint  satisfaction  problem  known 
as  “exact  cover  3,”  also  called  l-in-3  sat),  found  that  the 
complexity  scaled  polynomially  with  size,  roughly  as  N2, 
which  caused  a  good  deal  of  excitement.  However,  this 
power  law  complexity  may  be  an  artifact  of  the  very  small 
sizes  studied,  so  it  is  of  great  interest  to  determine  whether 
the  complexity  continues  to  be  polynomial  for  much  larger 
sizes  or  whether  a  “crossover”  to  exponential  complexity 
is  seen. 

In  previous  work  [5],  we  have  used  quantum 
Monte  Carlo  (QMC)  simulations  to  investigate  much  larger 
sizes  of  the  exact  cover  problem,  up  to  N  =  128.  We  found 
evidence  that,  while  the  median  complexity  is  still  poly¬ 
nomial,  an  increasing  fraction  of  instances  became  very 
hard  to  equilibrate  for  the  larger  sizes.  We  have  now 
considerably  improved  the  algorithm,  borrowing  tech¬ 
niques  from  the  spin  glass  field  to  speed  up  equilibration. 
We  have  therefore  been  able  to  understand  much  better 
these  “troublesome”  instances,  and  find  that  they  have  a 
first-order  quantum  phase  transition.  Furthermore,  we  have 
increased  the  range  of  sizes  still  further,  up  to  N  =  256, 
finding  that  the  fraction  instances  with  a  first-order  tran¬ 
sition  continues  to  increase  with  N,  plausibly  tending  1  for 
A  — ►  oo.  The  gap  at  a  first-order  phase  transition  is  likely  to 
be  exponentially  small  [6,7],  and  hence  lead  to  exponential 
complexity  for  the  QAA. 

We  now  describe  the  model  and  our  results  in  more 
detail.  To  make  a  comparison  with  the  earlier  work,  we 
study  (essentially)  the  same  model  for  3~C P  used  by  Farhi 
et  al,  [2].  This  problem,  known  as  exact  cover,  is  a  random 
satisfiability  problem,  a  class  which  is  known  to  be  NP 
hard.  In  exact  cover,  there  are  N  Ising  spins  and  M 
“clauses”  each  of  which  involves  three  spins  (chosen  at 
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random).  The  energy  of  a  clause  is  zero  if  one  spin  is  —  1 
and  the  other  two  are  1,  and  otherwise  the  energy  is  a 
positive  integer.  The  simplest  Hamiltonian  with  this  prop¬ 
erty  is  [8] 


J{P 


|  M 

t  X  +  +  ~  ^2’ 

4  a=  1 


(2) 


where  ax,  a2,  and  a 3  are  the  three  spins  in  clause  a,  and 

the  {&]}{i  =  I . A)  are  Pauli  matrices.  In  the  absence  of 

the  driver  Hamiltonian,  the  Pauli  matrices  can  be  replaced 
by  classical  Ising  spins  5(-  taking  values  ±  1 .  An  instance 
has  a  “satisfying  assignment”  if  there  is  at  least  one  choice 
for  the  spins  where  the  total  energy  is  zero.  As  the  ratio 
a  =  M/N  is  increased,  there  is  a  phase  transition  at  as 
where  the  number  of  satisfying  assignments  goes  to  zero. 
The  version  used  by  Farhi  et  al.  considers  only  instances 
with  a  unique  satisfying  assignment  (USA),  i.e.,  there  is 
only  one  state  with  energy  0.  This  has  the  advantage  that 
the  gap  A E(s)  between  the  ground  state  and  first  excited 
state  is  greater  than  zero  in  both  limiting  cases,  J~C  =  J~C D 
and  3~C  =  J~C P,  but  will  have  a  minimum  at  an  intermedi¬ 
ate  value  s  —  s* .  In  addition,  it  ensures  that  we  work  close 
to  the  satisfiability  transition  where  the  problem  is  particu¬ 
larly  hard  [9],  Hence  here,  and  in  the  earlier  work  [5],  we 
consider  instances  with  a  USA. 

The  method  of  generating  instances  with  a  USA  is 
described  in  Ref.  [5].  For  each  size  N,  we  choose  the 
number  of  clauses  M  which  maximizes  the  probability  of 
finding  a  USA,  see  Table  I.  The  actual  number  of  spins 
simulated  N',  is  somewhat  less  than  N  due  to  isolated  sites 
being  omitted,  and  others  that  do  not  affect  the  complexity 
are  also  “pruned  off”  [5].  The  value  of  a  =  M/N  seems  to 
be  close  to  the  critical  value  as  —  0.626  [10]  for  N  — ►  oo. 

In  QMC,  we  simulate  an  effective  classical  model  with 
Ising  spins  S,(r)  =  ±1  in  which  r  (0  <  t  <  (3  =  T~l)  is 
imaginary  time.  Following  common  practice,  we  discretize 
imaginary  time  into  LT  “time  slices”  each  representing 
At  =  /3/Lt  of  imaginary  time.  We  take  At  =  1. 

As  discussed  previously  [5],  for  (3AE »  1  (where 
A E  =  £j  —  E0  is  the  energy  gap),  and  r  <SC  /3,  the  time- 
dependent  correlation  function 

l  n'  lt 

C(r)  =  ^  X  X  +  t)5,(t0)),  (3) 

T  i=  1  'o=l 

is  a  sum  of  exponentials,  i.e., 

C(t)  =  q  +  XA«  exp[  — (£„  -  E0)t],  (4) 

n>  1 


TABLE  I.  The  sizes  studied  in  the  simulation. 


N 

16 

32 

64 

128 

192 

256 

M 

12 

23 

44 

86 

126 

166 

a 

0.7500 

0.7188 

0.6875 

0.6719 

0.6563 

0.6484 

where  the  An  are  constants  and  q,  the  spin  glass  order 
parameter,  is  given  by 


<7  =  V7  X<5'->2-  <5) 

i=l 

At  large  t,  the  sum  in  Eq.  (4)  is  dominated  by  the  term 
corresponding  to  the  first  excited  state  (n  =  1),  and  so 
A E  =  £j  —  E0  can  be  obtained  by  fitting  log[C(r)  —  g] 
against  t  for  large  t. 

We  have  considerably  improved  the  algorithm  relative  to 
that  in  Ref.  [5]  by  incorporating  “parallel  tempering” 
[11,12],  which  has  been  very  successful  in  speeding  up 
simulations  of  spin  glass  systems.  Whereas  in  spin  glasses, 
one  simulates  copies  of  the  system  at  different,  close-by 
temperatures,  in  the  quantum  case,  the  copies  are  at  differ¬ 
ent  values  of  the  control  parameter  s. 

As  already  mentioned,  the  focus  of  the  present  study  is 
to  determine  which  instances  have  a  first-order  transition. 
Parallel  tempering  is  very  good  at  equilibrating  the  system 
on  either  side  of  the  transition.  However,  it  is  still  difficult 
(i)  to  determine  exactly  where  the  transition  occurs,  be¬ 
cause  both  phases  are  metastable  in  the  region  where  they 
are  not  the  equilibrium  state,  and  (ii)  to  accurately  deter¬ 
mine  the  minimum  gap  for  first-order  instances,  because  it 
is  so  small.  We  have  performed  runs  starting  the  spins  both 
from  a  random  initial  configuration  and  from  the  solution 
of  the  problem  Hamiltonian.  If  we  start  by  “seeding”  the 
spins  with  the  exact  solution,  we  are  confident  that  the 
Monte  Carlo  simulation  is  in  the  correct  phase  for  5  close  to 
1.  It  is  also  in  the  correct  phase  for  small  5  because 
equilibration  is  easy  in  this  region.  Hence,  if  a  long  simu¬ 
lation  starting  the  spins  from  the  exact  solution  produces  a 
sharp  discontinuity,  we  feel  that  this  is  almost  certainly  the 
correct  behavior. 

In  order  to  investigate  whether  or  not  a  first-order  tran¬ 
sition  occurs,  we  focus  on  the  spin  glass  order  parameter  q 
defined  in  Eq.  (5).  The  expectation  value  of  q  is  always 
nonzero  because  of  terms  linear  in  the  &z  (magnetic  field 
terms)  in  the  Hamiltonian,  Eq.  (2).  To  determine  the  square 
of  the  average  without  bias,  we  simulated  two  copies  of  the 
spins  at  each  value  of  s  and  evaluate  (S/)2  as  {S //^{S A 
representative  result  for  an  instance  with  a  first-order  tran¬ 
sition  is  shown  in  Fig.  1 .  Note  the  very  rapid  increase  in  q 
over  a  very  small  range  of  s,  and  that  the  two  curves  on 
each  side  of  the  jump  are  obviously  displaced  vertically 
with  respect  to  each  other. 

The  dip  before  the  jump  clearly  seen  in  Fig.  1  provides 
clear  evidence  for  a  two-phase  coexistence,  and  hence  a 
first-order  transition,  for  the  following  reason.  If  both 
copies  are  in  the  same  phase,  then  the  mean  value  of  (5,-) 
is  the  same  in  both  copies.  However,  right  at  the  first-order 
transition,  one  copy  can  be  in  one  phase  (the  low-g  phase, 
say)  and  the  other  copy  in  the  other  (high-g)  phase.  The 
average  value  of  (5,-)  can  have  different  signs  in  the  two 
phases  for  some  sites  i.  Hence,  the  typical  Hamming 
distance  between  the  spin  configurations  in  the  two  copies 


020502-2 


PRL  104,  020502  (2010) 

0.6  — i — i — i — i — 


PHYSICAL  REVIEW  LETTERS 


week  ending 
15  JANUARY  2010 


FIG.  1  (color  online).  An  instance  with  a  first-order  transition 
with  N  =  128.  Note  the  expanded  horizontal  scale. 


can  be  even  greater  (and  so  q  even  smaller)  than  when  both 
copies  are  in  the  low- <7  phase.  In  every  instance  where  we 
observed  a  sharp  jump,  this  was  proceeded  by  a  dip.  Hence, 
we  use  the  dip  as  a  precise  criterion  for  a  first-order 
transition. 

Of  course,  even  a  first-order  transition  is  rounded  out  for 
a  finite-size  system.  To  estimate  the  size  of  the  rounding, 
we  need  to  consider  the  two  cases  A£min  »  T  and 
A£min  <5C  T  separately,  where  AEmin  is  the  minimum  value 
of  the  gap  at  the  transition.  If  A£jnin  »  T,  Ss  is  the  range 
of  .v  over  which  A  E  changes  by  an  amount  ALmm,  whereas 
if  AEjnin  <C  T,  Ss  is  the  range  of  s  over  which  A E  changes 
by  an  amount  equal  to  T.  Hence, 


f  AEmin^r)-1,  (A£mln  »  n 
m^r1,  (A Emin  «  T) 


Figure  2  shows  the  finite-size  rounding  for  an  instance  with 
N  =  64,  small  enough  that  we  can  equilibrate  through  the 
(first-order)  transition.  For  /3  £  1024,  the  width  of  the 
transition  region  increases  as  /3  =  1  /T  decreases,  but  for 
/ 3  S  1024  the  width  is  independent  of  [i.  For  this  instance, 
we  find  A  E^m  =  0.0021  as  shown  in  the  inset,  so  the  width 
of  the  rounding  becomes  independent  of  T  when  T  «  A E 
as  expected. 

In  Fig.  3,  we  plot  the  fraction  of  instances  with  a  first- 
order  transition.  For  each  size,  we  have  studied  V;nst  =  50 
instances.  If  we  denote  the  first-order  fraction  by  r,  then  the 
error  bar  in  r  is  *Jr(  I  —  r)/(Vjnst  —  1).  The  figure  shows 
that  r  increases  rapidly  with  N  and,  very  plausibly,  tends  to 
1  for  N  —*  oo.  We  see  that  the  first-order  fraction  is  slightly 
greater  than  a  half  for  N  =  128.  In  our  earlier  work  [5],  we 
found  that  the  median  complexity  continued  to  be  poly¬ 
nomial  up  to  N  =  128  (the  largest  size  studied).  However, 
there  is  no  contrast  with  the  present  work  because,  as 
already  noted  in  Ref.  [8],  the  models  used  are  slightly 


FIG.  2  (color  online).  The  main  figure  shows  the  spin  glass 
order  parameter  q,  defined  in  Eq.  (5),  as  a  function  of  s  for  an 
instance  with  N  =  64  which  has  a  first-order  transition.  The 
different  curves  are  for  different  values  of  p.  The  inset  shows  the 
energy  gap  Aft  as  a  function  of  s  for  /3  =  2048,  indicating  that 
Aftmin  =  0.0021  (same  value  was  found  for  /?  =  1024  and 
4096).  From  the  main  figure,  one  sees  that  the  width  of  the 
finite-size  rounding  increases  with  T  =  1/yS  for  T  »  Aft  but  is 
independent  of  T  in  the  opposite  limit  T  <3C  Aft,  as  expected 
from  Eq.  (6).  Note  the  expanded  horizontal  scale. 


different,  and  as  a  result,  the  crossover  to  a  first-order 
transition  occurs  at  a  slightly  lower  value  of  N  in  the 
present  model.  The  crossover  to  first  order  would  have 
been  seen  in  the  earlier  model  if  somewhat  larger  sizes 
had  been  studied. 

Exponentially  small  gaps  have  been  discussed  before  in 
the  context  of  the  QAA.  Some  time  ago,  one  of  us  [13] 


FIG.  3  (color  online).  The  fraction  of  instances  with  a  first- 
order  transition  (defined  in  the  way  discussed  in  the  text)  as  a 
function  of  size.  For  each  size,  50  instances  were  studied. 
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pointed  out  for  a  different  problem,  number  partitioning, 
that  the  minimum  gap  is  exponentially  small,  because  of  a 
transition  between  the  states  that  are  “localized”  and 
“extended”  in  the  computational  basis. 

Altshuler  et  al.  [6]  predict  an  exponentially  small  gap  at 
large  N  for  exact  cover.  Performing  perturbation  theory 
away  from  s  =  1,  they  argue  that  there  will  be  a  level 
crossing  between  two  “localized”  states  for  5  close  to  1, 
at  which  point  the  ground  state  configuration  changes 
abruptly.  In  our  numerics,  there  is  a  big  variation  in  the 
location  of  the  first-order  transition  for  a  given  size,  but  we 
do  not  detect  a  systematic  shift  towards  5  =  1  as  the  size 
increases.  However,  Altshuler  et  al.  predict  that  1  —  s  ~ 
V-1/8,  which  is  probably  too  slow  to  be  visible  in  our  data. 
It  will  be  interesting  to  investigate  in  future  work  whether 
the  first-order  transition  found  here  is  due  to  the  mecha¬ 
nism  they  propose. 

Farhi  et  al.  [14]  used  a  continuous  imaginary  time  QMC 
method  to  study  a  very  similar  problem  to  ours,  except  that 
two  solutions  far  away  in  Hamming  space  are  “planted” 
into  the  Hamiltonian.  This  ensures  that  there  is  a  finite 
probability  of  a  first-order  transition  where  the  equilibrium 
state  changes  from  one  planted  solution  to  another.  By 
contrast,  our  work  does  nothing  explicit  to  impose  a  first- 
order  transition. 

Jorg  et  al.  [15]  studied  quantum  annealing  for  the  quan¬ 
tum  random  energy  model  (REM),  the  classical  version  of 
which  [16]  has  a  “1-step  replica  symmetry  breaking”  (also 
called  a  “random  first-order”)  transition.  Following 
Goldschmidt  [17],  they  find  a  discontinuous  quantum  tran¬ 
sition  and  argue  that  this  leads  to  an  exponentially  small 
gap.  They  also  observed  that  an  exponentially  small  gap  is 
seen  in  quantum  versions  of  several  models  with  random 
first-order  transitions  and  suggested  that  this  may  be  the 
general  feature  of  all  such  models,  including  satisfiability 
[10,18],  However,  the  classical  REM  has  zero  spin  glass 
order  parameter  q  in  the  disordered  phase  [16]  whereas 
classical  random  satisfiability  models  have  lower  symme¬ 
try  because  q  is  always  nonzero  due  to  the  terms  linear  in 
&z  in  Eq.  (2).  Consequently,  it  is  not  obvious  to  us  that  the 
first-order  quantum  transition  observed  here  is  due  to  the 
same  mechanism  as  that  found  [15]  for  the  quantum  REM. 
Very  recently,  a  first-order  transition  has  also  been  found  in 
another  model  by  Jorg  et  al.  [19]. 

To  conclude,  we  have  a  found  a  crossover  to  a  first-order 
quantum  phase  transition  during  the  evolution  of  the  QAA 
for  instances  of  exact  cover  with  a  unique  satisfying  assign¬ 
ment  when  the  size  becomes  greater  than  about  100.  It  is 
possible  that  the  complexity  for  random  instances  of  exact 
cover  could  be  different.  We  are  therefore  studying  instan¬ 
ces  of  exact  cover  with  the  USA  constraint  removed,  and 
will  also  study  other  models  in  addition  to  exact  cover. 
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